#include "./cells-dyn.h"

#if defined(_OPENMP)
#include <omp.h>
#endif

void mds_cells_dynamics_whole(mds_frame         *frame, 
							  mds_potential_proc Vproc, 
							  void              *Vargs,
							  void              *params)
{
	
	assert(frame!=NULL);
	assert(frame->num_particles>0);
	assert(frame->num_sets>=MDS_ACCELERATION);
	assert(Vproc!=NULL);
	assert(params!=NULL);
	{
		mds_vector_set     positions    = frame->sets[MDS_POSITION];
		mds_vector_set     forces       = frame->sets[MDS_ACCELERATION];
		const mds_vector  *box          = & frame->box;
		mds_pbc_proc       pbc_for_f    = frame->pbc_for_f;
		mds_cells         *cells        = params;
		mds_pnode         *nodes        = frame->nodes;
		mds_cell          *origin       = & cells->array[0][0][0];
		const size_t       nCells       = cells->count;
		
		assert( cells != NULL );
		assert( nodes != NULL );
		
		/***************************/
		/** Initialize forces    **/
		/*************************/
		mds_frame_no_forces(frame);
		
		
		/****************************/
		/** dispatch particles    **/
		/**************************/
		mds_cells_dispatch( cells, positions, nodes, frame->num_particles );
		
		/**********************************************************************/
		/** outer loops on cells                                             **/
		/**********************************************************************/
		#pragma omp parallel default(shared) 
		{
#if defined(_OPENMP)
			const size_t __tid = omp_get_thread_num();
#else
#define                  __tid 0
#endif
			size_t       c;
			const size_t iEnter = cells->enter[ __tid ];
			const size_t iLeave = cells->leave[ __tid ];
			for( c=iEnter; c < iLeave; ++c )
			{
				Real V = 0;
				Real W = 0;
				Real V_ij, W_ij;
				mds_pnode         *node_i;
				mds_pnode         *node_j;
				mds_cell          *qCell = NULL;
				mds_cell          *pCell = &origin[c];
				mds_vector         f_ij;
				mds_vector         r_ij;
				mds_cell         **neighbors  = pCell->neighbors;
				const size_t       nq         = pCell->num_neighbors;
				size_t             q;
				
				
				for( node_i = pCell->first; node_i != NULL; node_i = node_i->next )
				{
					const size_t i   = node_i->id;
					mds_vector  *r_i = &positions[i];
					mds_vector  *f_i = &forces[i];
					
					/*------------------------------------------------------------*/
					/* Inside the cell                                            */
					/*------------------------------------------------------------*/
					for( node_j = node_i->next; node_j != NULL; node_j = node_j->next )
					{
						const size_t j   = node_j->id;
						mds_vector  *r_j = &positions[j];
						mds_vector  *f_j = &forces[j];
						assert( j != i );
						
						/** compute pair separation */
						mds_vector_vec( &r_ij, r_i, r_j );
						
						/** apply PBC for pair separation */
						pbc_for_f( &r_ij, box );
						
						/** update V and compute f_ij */
						V_ij = Vproc( &r_ij, &f_ij, Vargs );
						V += (V_ij + V_ij );
						
						/** dispatch forces */
						mds_vector_sub( f_i, &f_ij );
						mds_vector_add( f_j, &f_ij );
						
						/** update Virial */
						W_ij = mds_vector_dot( &r_ij, &f_ij );
						W += (W_ij + W_ij);
					}
					
					
					/*------------------------------------------------------------*/
					/* with the neighbors                                         */
					/*------------------------------------------------------------*/
					for( q=0; q < nq; ++q )
					{
						qCell = neighbors[q]; assert( qCell != NULL );
						for( node_j = qCell->first; node_j != NULL; node_j = node_j->next )
						{
							const size_t j   = node_j->id;
							mds_vector  *r_j = &positions[j];
							
							assert( i != j );
							
							/** compute pair separation */
							mds_vector_vec( &r_ij, r_i, r_j );
							
							/** apply PBC for pair separation */
							pbc_for_f( &r_ij, box );
							
							/** update V and compute f_ij */
							V += Vproc( &r_ij, &f_ij, Vargs );
							
							/** dispatch only one force ! */
							mds_vector_sub( f_i, &f_ij );
							
							/** update Virial */
							W += mds_vector_dot( &r_ij, &f_ij );						
						}
					}
					
				}
				
				pCell->V = V;
				pCell->W = W;
				
			}
		}
		
		/***********************************************************************
		 *   partial reduction             
		 **********************************************************************/
		{
			Real      V     = 0;
			Real      W     = 0;
			mds_cell *pCell = origin;
			register size_t c;
			for(c=0;c<nCells;++c,++pCell)
			{
				V    += pCell->V;
				W    += pCell->W;
			}
			
			frame->V = MDS_F(0.5) * V;
			frame->W = (MDS_F(0.5) * W) / MDS_F(3.0);
		}
		
	}
}
